Chapter 13 

Introduction to Simulation Using R 

A. Rakhshan and H. Pishro-Nik 


13.1 Analysis versus Computer Simulation 

A computer simulation is a computer program which attempts to represent the real world based 
on a model. The accuracy of the simulation depends on the precision of the model. Suppose that 
the probability of heads in a coin toss experiment is unknown. We can perform the experiment 
of tossing the coin n times repetitively to approximate the probability of heads. 

prjj\ Number of times heads observed 

Number of times the experiment executed 

However, for many practical problems it is not possible to determine the probabilities by exe¬ 
cuting experiments a large number of times. With today’s computers processing capabilities, we 
only need a high-level language, such as R, which can generate random numbers, to deal with 
these problems. 

In this chapter, we present basic methods of generating random variables and simulate 
probabilistic systems. The provided algorithms are general and can be implemented in any 
computer language. However, to have concrete examples, we provide the actual codes in R. If 
you are unfamiliar with R, you should still be able to understand the algorithms. 

13.2 Introduction: What is R? 

R is a programming language that helps engineers and scientists find solutions for given statisti¬ 
cal problems with fewer lines of codes than traditional programming languages, such as C/C-l—I- 
or Java, by utilizing built-in statistical functions. There are many built-in statistical functions 
and add-on packages available in R. It also has high quality customizable graphics capabilities. 
R is available for Unix/Linux, Windows, and Mac. Besides all these features, R is free! 

13.3 Discrete and Continuous Random Number Generators 

Most of the programming languages can deliver samples from the uniform distribution to us 
(In reality, the given values are pseudo-random instead of being completely random.) The rest 
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of this section shows how to convert uniform random variables to any other desired random 
variable. The R code for generating uniform random variables is: 


U = runif(n, min = 0, max = 1) 


which returns a pseudorandom value drawn from the standard uniform distribution on the open 
interval (0,1). 


13.3.1 Generating Discrete Probability Distributions from Uniform Distri¬ 
bution 

Let’s see a few examples of generating certain simple distributions: 


Example 1. (Bernoulli) Simulate tossing a coin with probability of heads p. 

Solution: Let U be a Uniform(0,l) random variable. We can write Bernoulli random variable 
X as: 


f 1 U < p 
\ 0 U >p 


Thus, 


P(H) = P(X = 1) 
= P(U < p) 

= p 


Therefore, X has Bernoulli(p ) distribution. 


The R code for Bernoulli (0.5) is: 


p = 0.5; 

U = runif(l,min = 0 ,max = 1); 
X = {U<p)- 


Since the “runif(l, min = 0, max = 1)” command returns a number between 0 and 1, we divided 
the interval [0,1] into two parts, p and 1 — p in length. Then, the value of X is determined based 
on where the number generated from uniform distribution fell. 

Example 2. (Coin Toss Simulation) Write codes to simulate tossing a fair coin to see how the 
law of large numbers works. 

Solution: You can write: 
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n = 1000; 

U = runif(n , min = 0, max = 1); 
toss = U < 0.5; 
a = numeric(n + 1); 
avg = numeric(n ); 
for(i in 2 : n + 1) 

{ 

a[i\ = a[i — 1] + toss[i — 1]; 
avg[i - 1] = a[i\/(i - 1); 

} 

plot( 1 : n,avg[l : n\,type = "V’,lwd = 5, col = ”blue”,ylab = ”Proportionsf Heads ”, 
xlab = ” CoinTossNumber” ,cex.main = 1.25, cex.lab = 1.5, cex.axis = 1.75) 


If you run the above codes to compute the proportion of ones in the variable “toss,” the result 
will look like Figure 13.1. You can also assume the coin is unbiased with probability of heads 
equal to 0.6 by replacing the third line of the previous code with: 


toss = U < 0.6; 



Figure 13.1: R - coin toss simualtion 
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Example 3. (Binomial) Generate a Binomial(50,0.2) random variable. 

Solution: To solve this problem, we can use the following lemma: 

Lemma 1. If X\,X 2 , • • • • X n are independent Bernoulli(p) random variables, then the random 
variable X defined by X = X\ + X 2 + ... + X n has a Binomial(n,p) distribution. 

To generate a random variable X ~ Binomial(n,p), we can toss a coin n times and count the 
number of heads. Counting the number of heads is exactly the same as finding X\ + X 2 +...+ X n , 
where each X, t is equal to one if the corresponding coin toss results in heads and zero otherwise. 

Since we know how to generate Bernoulli random variables, we can generate a Binomial(n,p) 
by adding n independent Bernoulli(p ) random variables. 


p = 0- 2 ; 
n = 50; 

U = runif(n, min = 0, max = 1); 
X = sum(U < p)\ 


Generating Arbitrary Discrete Distributions 


In general, we can generate any discrete random variables similar to the above examples using 
the following algorithm. Suppose we would like to simulate the discrete random variable X with 
range Rx = {x\,x 2 , ... ,x n } and P(X = Xj) = pj, so YljPj = 1 - 

To achieve this, first we generate a random number U (i.e., U ~ Uniform( 0,1)). Next, we 
divide the interval [0,1] into subintervals such that the jth subinterval has length pj (Figure 
13.2). Assume 


X= { 


Xo 

if 

(u < Po) 

X\ 

if 

(po<U < po+pi) 

Xj 

if 

(Efe=0 Pk < u < E 


In other words 

X = Xj if F(xj- 1) <U< F(xj ), 
where F(x) is the desired CDF. We have 


P(X = Xj ) = P 


3 -1 3 \ 

^2pk < u < 

,k =0 k =0 / 


= pj 
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|P0^1 | P2 | P3 | ••• { Pj | _ 

0 1 

Figure 13.2: Generating discrete random variables 


Example 4. Give an algorithm to simulate the value of a random variable X such that 

P(X = 1) = 0.35 
P(X = 2) = 0.15 
P(X = 3) = 0.4 
P(X = 4) = 0.1 

Solution: We divide the interval [0,1] into subintervals as follows: 

A 0 = [0,0.35) 

Ai = [0.35,0.5) 

A 2 = [0.5, 0.9) 

A 3 = [0.9,1) 

Subinterval Aj has length p^. We obtain a uniform number U. If U belongs to A;, then X = Xj. 

P(X = Xi ) = P{U g Ai) 

= Pi 


P = c(0.35, 0.5, 0.9,1); 

X = c(l, 2,3,4); 
counter = 1; 

r = runif( 1, min = 0, max = 1); 
while(r > P[counter \) 
counter = counter + 1; 
end 

X [counter\ 


13.3.2 Generating Continuous Probability Distributions from the Uniform 
Distribution- Inverse Transformation Method 

At least in principle, there is a way to convert a uniform distribution to any other distribution. 
Let’s see how we can do this. Let U ~ Uniform(0,l) and F be a CDF. Also, assume F is 
continuous and strictly increasing as a function. 
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Theorem 1. Let U ~ Uniform( 0,1) and F be a CDF which is strictly increasing. Also, 
consider a random variable X defined as 


X = F~ 1 (U). 


Then, 


Proof: 


A ~ F (The CDF of X is F) 


P(X < x ) = P{F^ 1 (U) < x) 

= P(U < F(x)) (increasing function) 
= F(x) 


Now, let’s see some examples. Note that to generate any continuous random variable X with 
the continuous cdf F, F _1 (C7) has to be computed. 

Example 5. (Exponential) Generate an Exponential(l) random variable. 

Solution: To generate an Exponential random variable with parameter A = 1, we proceed 
as follows 


F(x) = 1 — e x x > 0 
U Uniform,(0, 1) 

X = F~ l (U) 

= - ln(l - U) 

X ~ F 


This formula can be simplified since 


1 — U ~ Uniform{ 0,1) 


U , 1 - U 


Figure 13.3: Symmetry of Uniform 


Hence we can simulate X using 


X = - In (U) 


U = runif(l,min = 0, max = 1); 
X = —log(U) 
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Example 6. (Gamma) Generate a Gamma(20,l) random variable. 

Solution: For this example, F -1 is even more complicated than the complicated gamma cdf 
F itself. Instead of inverting the CDF, we generate a gamma random variable as a sum of n 
independent exponential variables. 

Theorem 2. Let X \, X 2 , • • • , X n be independent random variables with AQ ~ Exponential A). 
Define 


Y = X 1 +X 2 + --- + X n 


By the moment generating function method, you can show that Y has a gamma distribution 
with parameters n and A, i.e., Y ~ Gamma(n , A). 

Having this theorem in mind, we can write: 


n = 20; 
lambda = 1; 

X = {—1/lambda) * sum{log{runif{n,min = 0 ,max = 1))); 


Example 7. (Poisson) Generate a Poisson random variable. Hint: In this example, use the fact 
that the number of events in the interval [0, t] has Poisson distribution when the elapsed times 
between the events are Exponential. 

Solution: We want to employ the definition of Poisson processes. Assume N represents the 
number of events (arrivals) in [0,t]. If the interarrival times are distributed exponentially (with 
parameter A) and independently, then the number of arrivals occurred in [0,t], N, has Poisson 
distribution with parameter A t (Figure 13.4). Therefore, to solve this problem, we can repeat 
generating Exponential{X) random variables while their sum is not larger than 1 (choosing 
t = 1). More specifically, we generate Exponential{X) random variables 

T % = ^ ln(?7j) 

by first generating uniform random variables U/'s. Then we define 


X = max { j '■ T\ + ■ ■ ■ + T.j < 1} 


The algorithm can be simplified: 




max 


-1 

X 


111 ( 1 /! • 


Uj)< 1 
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Lambda = 2; 
i = 0; 

U = runif(l,min = 0 ,max = 1); 

Y = —{1/Lambda) * log(U ); 
sum = Y ; 

while(sum < 1) 

{U = runif( 1, min = 0, max = 1); 

Y = —{1/Lambda) * log(U ); 
sum = sum + Y ; 

i = i + 1; } 

X = * 


I Exp(X) ! Exp(X) j Exp(X) ! Exp(X) ^ | 

0 1 
X= Maximum number of exponential random variables 

Figure 13.4: Poisson Random Variable 


To finish this section, let’s see how to convert uniform numbers to normal random variables. 
Normal distribution is extremely important in science because it is very commonly occuring. 


Theorem 3. (Box-Muller transformation) We can generate a pair of independent normal vari¬ 
ables (Z i, Z- 2 ) by transforming a pair of independent Uniform(Q , 1) random variables {U\, U 2 ) 
[!]• 


J Z\ = yj — 2 In U\ COs(27t[/2) 
I Z 2 = \/—2 In U\ sin(27rt/2) 


Example 8. (Box-Muller) Generate 5000 pairs of normal random variables and plot both 
histograms. 

Solution: We display the pairs in Matrix form. 
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n = 5000; 

U 1 = runif(n,min = 0 ,max = 1) 

U 2 = runif(n, min = 0, max = 1) 

Z 1 = sqrt(— 2 * log(Ul )) * cos(2 * pi * U2) 
Z2 = sqrt(—2 * log{U \)) * sin(2 * pi * U2) 
hist(Zl,col = ” wheat ", label = T) 


Histogram of Z1 



Z1 


Figure 13.5: Histogram of Zl, a normal random variable generated by Box-Muller transforma¬ 
tion 


13.4 R Commands for Special Distributions 

In this section, we will see some useful commands for commonly employed distributions. Func¬ 
tions which start with “p”,“q”,“d”, and “r” give the cumulative distribution function (CDF), 
the inverse CDF, the density function (PDF), and a random variable having the specified dis¬ 
tribution respectively. 


Distributions 


Commands 
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Binomial 

pbinom. 

qbinom 

dbinom 

rbinom 

Geometric 

pgeom 

qgeom. 

dgeom 

rgeom. 

NegativeBinomial 

pnbinom 

qnbinom 

dnbinom 

rnbinom 

Poisson 

ppois 

qpois 

dpois 

rpois 

Beta 

pbeta 

qbet.a 

dbeta 

rbeta 

Beta 

pbeta 

qbeta 

dbeta 

rbeta 

Exponential 

pexp 

qexp 

dexp 

rexp 

Gamma 

pgamma 

qgamma 

dgamma 

rgamma 

Studentt 

pt 

qt 

dt 

rt 

Uniform 

punif 

qunif 

dunif 

runif 


13.5 Exercises 


1. Write R programs to generate Geometric(p) and Negative Binomial(i,p) random variables. 

Solution: To generate a Geometric random variable, we run a loop of Bernoulli trials until 
the first success occurs. K counts the number of failures plus one success, which is equal 
to the total number of trials. 


K= 1; 

p = 0.2; 

while(runif( 1) > p) 
K = K + 1; 

K 


Now, we can generate Geometric random variable i times to obtain a Negative Binomial (i,p) 
variable as a sum of i independent Geometric (p) random variables. 
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2. (Poisson) Use the algorithm for generating discrete random variables to obtain a Poisson 
random variable with parameter A = 2. 

Solution: We know a Poisson random variable takes all nonnegative integer values with 
probabilities 

Pi = P(X = Xi) = e ~ x -- for i = 0,1,2, ••• 

i\ 

To generate a Poisson{ A), first we generate a random number U. Next, we divide the 
interval [0,1] into subintervals such that the jth subinterval has length pj (Figure 13.2). 
Assume 

x 0 if (U < po ) 

xi if (jpq < U < po+Pi) 

X = 

x j if (Ei=l'Pk <u < Yi=oPk) 


Here Xi = i — 1, so 


X = i if po H - Pi-i < U < po -\ -h pi-i + Pi 

F(i - 1) < U < F{i) F is CDF 
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3. Explain how to generate a random variable with the density 

f(x) = 2.5xy / x for 0 < x < 1 

if your random number generator produces a Standard Uniform random variable U. Hint: 
use the inverse transformation method. 

Solution: 


F x (X) = xl=U (0 < x < 1) 

X = ui 



We have the desired distribution. 
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By generating a random number, U, we have the desired distribution. 


U = runif{ 1); 
x = S( l rt ( 2U+ l) 


5. Let X have a standard Cauchy distribution. 

Fx (x) = — arctan(x') + - 
7r 2 


Assuming you have U ~ Uniform^ 0,1), explain how to generate X. Then, use this result 
to produce 1000 samples of X and compute the sample mean. Repeat the experiment 100 
times. What do you observe and why? 

Solution: Using Inverse Transformation Method: 

U -= — arc tan (A) 

2 7T 

7T (u — = arctan(A) 

X = tan ^7 t(U — -)\ 


Next, here is the R code: 


U = numeric( 1000); 
n = 100; 

average = numeric{n)\ 
for (i in 1 : n) 

{U = runif( 1000); 

X = tan(pi * (U — 0.5)); 
average[i\ = mean(X);} 

plot (l : n, average[l : n\,type = '"Vflwd = 2 ,col = "blue") 


Cauchy distribution has no mean (Figure 13.6), or higher moments defined. 

6. (The Rejection Method) When we use the Inverse Transformation Method, we need a 
simple form of the cdf F(x) that allows direct computation of X = I ?_1 (C/). When 
F(x) doesn’t have a simple form but the pdf f(x) is available, random variables with 
density f(x) can be generated by the rejection method. Suppose you have a method 
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Figure 13.6: Cauchy Mean 


for generating a random variable having density function g(x). Now, assume you want to 
generate a random variable having density function f{x). Let c be a constant such that 


< c (for all y) 

a{y) ~ 


Show that the following method generates a random variable with density function f(x). 

- Generate Y having density g. 

- Generate a random number U from Uniform (0,1). 

- If U < cg(y) , set X = Y. Otherwise, return to step 1. 

Solution: The number of times N that the first two steps of the algorithm need to be called 
is itself a random variable and has a geometric distribution with “success” probability 


p = P 



f(Y) \ 

cg(Y)J 
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Thus, E(N) = Also, we can compute p: 


P U < 


f(X) 
cg(Y ) 


Y = y\ = 


f(y ) 

cg(y ) 


r m , w 

P = / —r^9{y)dy 

7-oo w) 

1 Z 100 

= - f(y)dy 

( ' J—oo 


Therefore, E(N) = c 

Let F be the desired CDF (CDF of A). Now, we must show that the conditional distri¬ 
bution of y given that U < ^^ is indeed F, i.e. P(y < y|J7 < ) = F(y). Assume 

M = {U < cg(Y) }> K = {Y < y}. We know P(M) = p = I. Also, we can compute 

P(U<m,Y<y) 


P{U < PY\Y < y) = 
cg{Y) 


Thus, 


_ cg(Y )' _ 

G(y) 

ry P(U < $$\Y = v < y) 


' —OO 

1 


G(y) 


g(v)dv 


W) J-oo^) 9iv)dv 

^c)L mdv 

F(y) 

cG(y ) 


P(K\M) = P(M\K)P(K)/P{M) 

=p(us !$j' 

F(y) v G(y) 


cG(y) 

= ^(y) 


7. Use the rejection method to generate a random variable having density function Beta(2,4). 
Hint: Assume g(x) = 1 for 0 < x < 1. 

Solution: 

/(x) = 20x(l — x ) 3 0 < x < 1 

g(x) = 1 0 < x < 1 

= 20x(l - x) 3 

g{x) 
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We need to find the smallest constant c such that f{x)/g(x ) < c. Differentiation of this 
quantity yields 


/ 0) 
flU) 


dx 


= 0 


Thus, x = - 
4 


Therefore, 


Hence, 


fix) < 135 
g{x ) 
fix) 


64 


Cfi((x) 


256 n ^ 

= -^i 1 -*) 


n = 1; 

while fn == 1) 

{C/1 = runifi 1); 

1/2 = runifi 1); 

if {U2 <= 256/27* HI * (1 - t/l) A 3) 

n = 0;}} 


8. Use the rejection method to generate a random variable having the Gamma {§, 1) density 
function. Hint: Assume g(x) is the pdf of the Gamma (a = |, A = l). 

Solution: 


fix) 

9ix) 



x > 0 


2 __2x 

-e s x > 0 
5 


/(g) _ 

9(g) 3^vr 


10 3 

x2 e 


d(/& 

SW 


dx 

Hence, x = 


= 0 

5 
2 


10 / 5 

C “ VT U 


3rr 

5 


3 

2 =3 

e 2 


/(g) 

cg(x) 


3 

X2e 



— 3x 
5 


-3 

e 2 


We know how to generate an Exponential random variable. 
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- Generate a random number U\ and set Y = — | \ogU\. 

- Generate a random number Ui- 


-If U 2 < 


3 -3 Y 

Y le 


- 3 J _ 3 , set X = Y. Otherwise, execute the step 1. 

(§) 5 ^ 


9. Use the rejection method to generate a standard normal random variable. Hint: Assume 
g(x ) is the pdf of the exponential distribution with A = 1. 

Solution: 


/(*) = 


2 

e 2 


\/27T 


0 < .t < oo 


Thus 


5 <x) = e 
f{x 


0 < x < oo (Exponential density function with mean 1) 


2 z-2 2 

= \l ~e 2 

7T 


g(x) 

Thus, x = 1 maximizes 

Thus, c = 

/(*) 


f(x) 

g{x) 


(£-ir 




= e 


Generate T, an exponential random variable with mean 1. 
Generate a random number U. 


- If 17 < e 


-(y-ir 


set A = Y. Otherwise, return to step 1. 


10. Use the rejection method to generate a Gamma(2, 1) random variable conditional on its 
value being greater than 5, that is 


f(x) = 


xe 


roc 

Ji 5 Xe 


xe 


6e -5 


~ x dx 
{x > 5) 


Hint: Assume g(x) be the density function of exponential distribution. 

Solution: Since Gamma(2,l) random variable has expected value 2, we use an exponential 
distribution with mean 2 that is conditioned to be greater than 5. 

qf.pl x ) 


f(x) 

- f 5 °°xe(~ 

' X ldx 


xe(~ x ) 

x > 5 


6e( -5 ) 

g(x) 

l e (-S) 

2 e 

-5 

x > 5 

/(*) 

e 2 

X (x — 5 

= — e v 2 

) 


g(x) 
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We obtain the maximum in x = 5 since is decreasing. Therefore, 

c= m = 5_ 

5(5) 3 

- Generate a random number V. 

- Y = 5-21og(G). 

- Generate a random number U. 

- If U < 2 \ set X = Y: otherwise return to step 1. 
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